      subroutine read_paramfile(paramfile,maxlev,ndim_fft,whichrand,
     &  kmin_out,kmax_out,sphmode,cic_correct,iquad,lbox,downweight,
     &  xcen,ycen,zcen,side,sigma8_p,zshift,omega0,lambda0,rnseed,nlevels,
     &  partfile,binary_ic_file,psfile)
      implicit none

      include 'mpif.h'
      integer IERR, MYID, NUMPROCS
      integer status(MPI_STATUS_SIZE)

      character*(*) paramfile
      integer*8 maxlev

      real twopi
      integer*8 itag
      parameter (itag=1003004009)
      parameter(twopi=2.*3.1415926535)     

      integer len_trim

      integer*8 iresim

      real*8    xcen(maxlev),ycen(maxlev),zcen(maxlev),side(0:maxlev)
      real*8 lbox(maxlev),downweight(maxlev)
      integer*8 idep(maxlev),ndim_fft(maxlev),nlevel(maxlev)
      integer*4 rnseed(maxlev)
      integer*8 whichrand(maxlev)
      integer*8 kmin_out(maxlev),kmax_out(maxlev),sphmode(maxlev)
      integer*8 cic_correct(maxlev),iquad(maxlev),nlevels

      real*8 llbox,omega0,lambda0,sigma8_p
      real*8 xc,yc,zc,sz
      real*8 zshift,aexp,sig8,lingro
      real*8 rkmin,scale

      integer*8 icheck,npart,ns
      integer*8 kmin,kmax,kmax_p
      integer*8 itag_a
      integer*8 iparent
      integer*8 i,j,iloop

      integer*4 nlev_4

      character*512   binary_ic_file
      character*512   partfile,psfile


      call MPI_COMM_SIZE(MPI_COMM_WORLD, NUMPROCS,IERR)
      call MPI_COMM_RANK(MPI_COMM_WORLD, MYID, IERR)

      if (MYID.eq.0) then !-------------------------------------Processor 0 only ---------

      open (10,file=paramfile,status='old')
      read (10,*) itag_a          ! Check for version compatibility
      if (itag.ne.itag_a) then
          stop 'Version number conflict. Parameter file likely to be inapplicable '
      endif
      read (10,*)!-----------------------------------------------------
      read (10,*) iresim
      if (iresim.ne.1) stop 'Only resimulations supported '
      read (10,*)!-----------------------------------------------------
      read (10,*)!*********** 
      read (10,*)!-----------------------------------------------------
      read (10,*)  partfile 
      read (10,*)  binary_ic_file ! Name of binary output file
      read (10,*)!-----------------------------------------------------
      read (10,*)  psfile
      read (10,*)!-----------------------------------------------------
      read (10,*)!************ Details of the parent simulation
      read (10,*) llbox           ! Parent box size in Mpc/h
      read (10,*) omega0,lambda0  ! Cosmological parameters of parent sim
      read (10,*) sigma8_p        ! sigma_8(z=0) of parent simulation
      read (10,*) zshift          ! Redshift of the initial conditions (bin file)

      do i=1,maxlev 
        read (10,*)!---------------------------------------------------
        read (10,*,end=42) icheck
        if (icheck.ne.i) stop 'Incorrect labelling of levels!'
        read (10,*) xcen(i)
        read (10,*) ycen(i)
        read (10,*) zcen(i)
        read (10,*) side(i)
        read (10,*) nlevel(i)       ! Number of particles in parent simulation
        read (10,*) rnseed(i)       ! Random number seed used in parent run.
        read (10,*) whichrand(i)    ! Which way of setting random phases
        read (10,*) ndim_fft(i)     ! Dimension of fourier mesh 
        read (10,*) cic_correct(i)  ! Correct disp to compensate assign sch
        read (10,*) iquad(i)        ! Quad interp - if so disable cic option
        read (10,*) idep(i)         ! Number of parent fourier grid
      enddo
 42   continue
      close(10)


      nlevels = i-1                  ! Number of fourier meshes   
      print*,'Number of fourier meshes = ',nlevels


      call check_nesting(nlevels,xcen,ycen,zcen,sphmode,side(1),idep)


      open (12,file='wave_setup.log',status='unknown')

   

      side(0) = 1.0   


c----------  Loop over fourier meshes  -----------


      do iloop = 1,nlevels

        lbox(iloop)  = side(iloop)*llbox 
        ns = rnseed(iloop)
        npart = nlevel(3)   ! <=== redundant
        xc = xcen(iloop)
        yc = ycen(iloop)
        zc = zcen(iloop)
        sz = side(iloop)
        iparent = idep(iloop)
        scale = side(iloop)/side(idep(iloop))

        if (iparent.eq.0) then
                  kmax_p = 1
        else
                  kmax_p = int((real(nlevel(iparent))+0.1)**(1./3.))/2 
        endif 

        kmax = int((real(nlevel(iloop))+0.1)**(1./3.))/2   ! Nyquist limit

        rkmin = max(0.0, scale*real(kmax_p-1))

        kmin = 1+rkmin
                 
        if (iparent.eq.0) kmin = 0

        kmin_out(iloop) = kmin
        kmax_out(iloop) = kmax


        print*,'Waves populated from ',kmin,  ' to (but not including)',kmax
        write (12,*) 'Waves populated from ',kmin,  ' to (but not including)',kmax

        if (iparent.eq.0) then
           downweight(iloop)  = 1.0
        else
           downweight(iloop)=sqrt(0.5+real(kmin)-scale*(real(kmax_p-1)+0.5))
        endif

        if (iquad(iloop).ne.0) cic_correct(iloop)=0   !  For quad interp make sure
                                                      !  no cic-correction is used.


        print*,iloop,' level waves '
        print*,'ns,lbox,npart',ns,lbox(iloop),npart
        print*,'ndim_fft ',ndim_fft(iloop)
        print*,'xc,yc,zc,sz ',xc,yc,zc,sz
        print*,'kmin,kmax ',kmin,kmax
        print*,'downweight ',downweight(iloop)
        print*

        write (12,*) iloop,' level waves '
        write (12,*) 'ns,lbox,npart',ns,lbox(iloop),npart
        write (12,*) 'ndim_fft ',ndim_fft(iloop)
        write (12,*) 'xc,yc,zc,sz ',xc,yc,zc,sz
        write (12,*) 'kmin,kmax ',kmin,kmax
        write (12,*) 'downweight ',downweight(iloop)
        write (12,*) 


      print*,'ns,lbox,npart,ndim_fft1 ',ns,lbox(iloop),npart,ndim_fft(iloop)
      print*,'xc,yc,zc,sz ',xc,yc,zc,sz
      print*,'kmin,kmax,downweight',kmin,kmax,downweight(iloop)




      enddo ! iloop


        
      close(12)  ! wave_direct.log


      endif  !   End processor-0 only.

c     Broadcast values to all processors


      call MPI_BCAST(nlevels,1,MPI_INTEGER8,0,MPI_COMM_WORLD,IERR)
      nlev_4 = nlevels    ! integer*4 parameter
      call MPI_BCAST(ndim_fft,nlev_4,MPI_INTEGER8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(whichrand,nlev_4,MPI_INTEGER8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(kmin_out,nlev_4,MPI_INTEGER8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(kmax_out,nlev_4,MPI_INTEGER8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(sphmode,nlev_4,MPI_INTEGER8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(cic_correct,nlev_4,MPI_INTEGER8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(iquad      ,nlev_4,MPI_INTEGER8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(lbox,nlev_4,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(downweight,nlev_4,MPI_REAL8,0,MPI_COMM_WORLD,IERR) 
      call MPI_BCAST(xcen,nlev_4,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(ycen,nlev_4,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(zcen,nlev_4,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(side(0),nlev_4+1,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(sigma8_p,1,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(zshift,1,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(omega0,1,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(lambda0,1,MPI_REAL8,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(rnseed,nlev_4,MPI_INTEGER,0,MPI_COMM_WORLD,IERR) 
      call MPI_BCAST(partfile,512,MPI_CHARACTER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(binary_ic_file,512,MPI_CHARACTER,0,MPI_COMM_WORLD,IERR)
      call MPI_BCAST(psfile,512,MPI_CHARACTER,0,MPI_COMM_WORLD,IERR)


 
c      do i=0,numprocs-1
c      if (MYID.eq.i) then
c             
c      print*,'MYID, nlevels = ',MYID,nlevels
c      print*,'MYID, ndim_fft =',MYID,ndim_fft(1)
c      print*,'MYID, whichrand =',MYID,whichrand(1)
c      print*,'MYID, kmin_out =',MYID,kmin_out(1)
c      print*,'MYID, kmax_out =',MYID,kmax_out(1)
c      print*,'MYID, sphmode =', MYID,sphmode(1)
c      print*,'MYID, cic_correct =',MYID,cic_correct(1)
c      print*,'MYID, lbox =',MYID,lbox(1)
c      print*,'MYID, downweight =',MYID,downweight(1)
c      print*,'MYID, xcen =',MYID,xcen(1)
c      print*,'MYID, ycen =',MYID,ycen(1)
c      print*,'MYID, zcen =',MYID,zcen(1)
c      print*,'MYID, side =',MYID,side(0),side(1)
c      print*,'MYID, sigma8_p',MYID,sigma8_p
c      print*,'MYID, rnseed',MYID,rnseed(1)
c      call flush(6)
c      endif
c            call MPI_BARRIER(MPI_COMM_WORLD,IERR)
c      enddo          
  
      

      end    





      subroutine check_nesting(numlev,xcen,ycen,zcen,sphmode,side,idep)
      implicit none
      integer*8 numlev
      real*8 xcen(numlev),ycen(numlev),zcen(numlev)
      real*8 side(numlev)
      real*8 dx,dy,dz,xsep,ysep,zsep
      integer*8 idep(numlev),sphmode(numlev)
      integer*8 icheck,iparent,icheck2
      integer*8 i

c  Default set sphmode = 1   -  for spherical boundary.  Set to zero for parent
c  of any level

      do i=1,numlev
         sphmode(i) = 1
      enddo

c Check top level

      icheck = 1

      if ((xcen(icheck).ne.0.5).or.(ycen(icheck).ne.0.5).or.(zcen(icheck).ne.0.5)) then
        print*,'Warning:  level 1 fourier box NOT CENTRED '
        print*,'xcen,ycen,zcen ',xcen(icheck),ycen(icheck),zcen(icheck)
        stop
      endif


      if (side(icheck).ne.1) then
           print*,'Warning: level 1 fourier box NOT of unit length'
           print*,'side(1) = ',side(icheck)
           stop
      endif 



      if (idep(icheck).ne.0) stop 'Level 1: has wrong parent'

c Now check subsequent levels



      do icheck=2,numlev
      print*,'icheck = ',icheck
c  Firstly parent box should be bigger than child box

        iparent = idep(icheck)
        sphmode(iparent) = 0     ! Cubic k-cut off for all parent boxes
        if (side(icheck).ge.side(iparent)) then
          print*,'level:',icheck,' fourier box greater or equal to parent!'
          print*, side(icheck),'>',side(iparent)
          print*,'icheck,iparent ',icheck,iparent
          stop
        endif


c  Box centre coordinates should be in principal region 0-1

        if ((xcen(icheck).lt.0).or.(xcen(icheck).ge.1)) then
           print*,'level:',icheck, ' xcen out of range 0-1 ',xcen(icheck)
           stop
        endif


        if ((ycen(icheck).lt.0).or.(ycen(icheck).ge.1)) then
           print*,'level:',icheck, ' ycen out of range 0-1 ',ycen(icheck)
           stop
        endif


        if ((zcen(icheck).lt.0).or.(zcen(icheck).ge.1)) then
           print*,'level:',icheck, ' zcen out of range 0-1 ',zcen(icheck)
           stop
        endif


c  Box of child should be within parent 
c    except if the parent is level 1 (and therefore periodic).

        if (iparent.ne.1) then


          dx = abs(xcen(icheck)-xcen(iparent))
          dy = abs(ycen(icheck)-ycen(iparent))
          dz = abs(zcen(icheck)-zcen(iparent))

          xsep = min(dx,1.-dx)
          ysep = min(dy,1.-dy)
          zsep = min(dz,1.-dz)


          if (0.5*side(iparent).lt.(xsep+0.5*side(icheck))) then
            print*,'level: ',icheck, ' child overlaps parent in x-direction '
            stop
          endif


          if (0.5*side(iparent).lt.(ysep+0.5*side(icheck))) then
            print*,'level: ',icheck, ' child overlaps parent in y-direction '
            stop
          endif


          if (0.5*side(iparent).lt.(zsep+0.5*side(icheck))) then
            print*,'level: ',icheck, ' child overlaps parent in z-direction '
            stop
          endif
        endif



      enddo
 

c   Check to see if boxes on the same level overlap 

                print*,'icheck, numlev ',icheck, numlev

      do icheck=1,numlev

       do icheck2=icheck+1,numlev

        if ((idep(icheck).eq.idep(icheck2))) then
           dx = abs(xcen(icheck)-xcen(icheck2))
           dy = abs(ycen(icheck)-ycen(icheck2))
           dz = abs(zcen(icheck)-zcen(icheck2))

           xsep = min(dx,1.-dx)
           ysep = min(dy,1.-dy)
           zsep = min(dz,1.-dz)



           if ((2.*xsep.le.(side(icheck)+side(icheck2))).and.
     &         (2.*ysep.le.(side(icheck)+side(icheck2))).and. 
     &         (2.*zsep.le.(side(icheck)+side(icheck2)))) then
             print*,'levels ',icheck,icheck2, ' overlap '
             stop
           endif
        endif
       enddo
      enddo    


      return
      end



